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Abstract.  We  review  time-domain  formulations  of  radiation  boundary  condi¬ 
tions  for  Maxwell’s  equations,  focusing  on  methods  which  can  deliver  arbitrary 
accuracy  at  acceptable  computational  cost.  Examples  include  fast  evaluations 
of  nonlocal  conditions  on  symmetric  and  general  boundaries,  methods  based  on 
identfying  and  evaluating  equivalent  sources,  and  local  approximations  such 
as  the  perfectly  matched  layer  and  sequences  of  local  boundary  conditions. 
Complexity  estimates  are  derived  to  assess  work  and  storage  requirements  as 
a  function  of  wavelength  and  simulation  time. 


1.  Introduction 

As  the  radiation  of  energy  to  the  far  field  is  an  important  feature  of  most  prob¬ 
lems  in  computational  electromagnetics,  an  accurate  and  efficient  truncation  of  the 
domain  is  a  practical  necessity  for  computations.  In  recent  years  there  have  been 
rapid  developments  in  this  field.  In  this  review  we  will  concentrate  on  strategies 
which  can  provide  arbitrary  accuracy.  These  include  a  variety  of  exact  boundary 
condition  formulations,  which  are  all  nonlocal  in  space  and  time,  in  addition  to  con¬ 
vergent  local  approximations  such  as  the  perfectly  matched  layer  (PML).  Besides 
describing  the  basic  mathematical  and  algorithmic  content  of  the  various  methods, 
we  will,  when  possible,  estimate  their  computational  complexity  as  a  function  of  the 
harmonic  content  of  the  field  and  the  simulation  time.  Our  goal  is  not  to  advocate 
one  of  the  methods  discussed  over  another.  We  will  see  that  they  are  all  capable 
of  providing  excellent  accuracy  at  acceptable  cost  in  many  settings,  and  that  an 
optimal  choice  will  depend  both  on  the  details  of  the  problem  as  well  as  on  the  time 
to  be  invested  on  code  development. 

We  will  assume  that  in  the  far  field,  that  is  beyond  the  computational  domain  f 1, 
we  have  a  homogeneous,  isotropic,  dielectric  material.  In  cgs  units  the  source-free 
Maxwell  equations  then  are: 


(1.1) 

dE 

— —  cVxB  =0, 
at 

(1.2) 

dB 

-j—  +  cV  x  E  —  0, 
at 

1991  Mathematics  Subject  Classification.  65M99,78M99. 

Supported  in  part  by  ARO  Grant  DAAD 19-03- 1-0146  and  NSF  Grant  DMS-0610067.  The 
second  author  was  also  supported  by  NSF  Grant  DMS-0554377.  Any  conclusions  or  recommen¬ 
dations  expressed  in  this  paper  are  those  of  the  authors  and  do  not  necessarily  reflect  the  views 
of  ARO  or  NSF. 


1 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

2007 


2.  REPORT  TYPE 


3.  DATES  COVERED 

00-00-2007  to  00-00-2007 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


4.  TITLE  AND  SUBTITLE 

Radiation  Boundary  Conditions  for  Maxwell’s  Equations:  A  Review  of 
Accurate  Time-Domain  Formulations 

6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  8.  PERFORMING  ORGANIZATION 

Brown  University, Division  of  Applied  Mathematics, 182  George  report  number 

Street, Providence, RI, 02912 

9.  SPONSORING/MONITORING  AGENCY  NAME(S )  AND  ADDRESS(ES )  10.  SPONSOR/MONITOR' S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

We  review  time-domain  formulations  of  radiation  boundary  conditions  for  Maxwell’s  equations,  focusing 
on  methods  which  can  deliver  arbitrary  accuracy  at  acceptable  computational  cost.  Examples  include  fast 
evaluations  of  nonlocal  conditions  on  symmetric  and  general  boundaries,  methods  based  on  identfying  and 
evaluating  equivalent  sources,  and  local  approximations  such  as  the  perfectly  matched  layer  and  sequences 
of  local  boundary  conditions.  Complexity  estimates  are  derived  to  assess  work  and  storage  requirements  as 
a  function  of  wavelength  and  simulation  time. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

33 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


2 


THOMAS  HAGSTROM  AND  STEPHEN  LAU 


subject  to  the  constraints 

(1.3)  V  •  E  =  0  =  V  •  B. 

The  constraints  (1.3)  are  clearly  preserved  under  the  time  evolution  governed  by 

(1-1.)  (I-2)- 

Our  problem  is  to  specify  radiation  boundary  conditions  at  an  artificial  boundary 
T  C  dS l  so  that  the  solution  computed  in  Q  can  be  made  arbitrarily  close  to  the 
restriction  to  17  of  the  solution  of  the  original  problem  on  the  unbounded  domain. 
We  will  organize  the  discussion  around  four  general  classes  of  methods:  fast  methods 
based  on  separation  of  variables  on  symmetric  boundaries,  methods  for  general 
boundaries  beased  on  the  retarded  potential,  methods  based  on  equivalent  source 
representations,  and,  finally,  convergent  local  approximations.  We  note  that  there 
have  been  parallel  developments  for  other  applications  and  refer  the  reader  to  [35, 
36]  for  more  comprehensive  if  slightly  older  reviews. 


2.  Boundaries  with  symmetry 

For  plane,  spherical,  and  cylindrical  boundaries,  this  section  formulates  exact 
nonreflecting  boundary  conditions  for  the  homogeneous  Maxwell  equations.  An 
earlier  review  article  [35]  described  these  boundary  conditions  from  a  more  general 
perspective.  In  contrast,  our  presentation  here  considers  only  the  Maxwell  system 
(1.1) — (1.2)  and  derives  the  relevant  boundary  conditions  and  effective  numerical 
approximations  from  the  ground  up. 


2.1.  Planar  boundary.  Let  X\  =  x  =  0  specify  the  planar  boundary  of  the  “com¬ 
putational  domain”  x  <  0.  On  the  system  (1.1) — (1.2)  we  perform  both  a  Laplace 
transform  (denoted  by  a  hat)  in  time  and  a  Fourier  transform  (denoted  by  a  bar)  in 
the  tangential  variables  (x2,x3)  =  ( y,z ),  thereby  obtaining  a  differential-algebraic 
system.  With  (k2,  £3)  representing  the  Fourier  variables  dual  to  (y,  z),  the  system’s 
algebraic  sector  is 


(2.4) 


sEi  =  ik2B3  —  ik3B2,  sB  1  =  ~ik2E3  +  ik3E2, 


where  s  =  s/c.  Using  these  algebraic  equations,  we  may  then  express  the  remaining 
differential  sector  solely  in  terms  of  the  tangential  variables  as  follows: 


(2.5) 
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The  eigenvalues  of  the  matrix  are 

(2.6)  A±  =  iry/.s2  •  14  -  14  =  ±  v/s2  +  |fc|2, 


with  each  one  doubly  degenerate.  In  (2.6)  we  define  the  branch  to  ensure  that 
A+  has  positive  real  part  for  Res  >  0,  and  A+  ~  s  as  s  — ■>  00.  As  the  branch 
cut  we  choose  a  curve  in  the  left-half  s- plane  running  from  i\k\  to  —  i\k\.  Our 
radiation  conditions  demand  that  solutions  to  (2.5)  remain  bounded  as  x  — »  00. 


RADIATION  CONDITIONS  FOR  MAXWELL’S  EQUATIONS 


3 


Such  solutions  are  then  of  the  form 
(2.7) 


[a(s,  k2,k3) 

(  s2  +  k2  \ 
k2k3 

0 

+  b(s,  k2,  k3) 

(  ~k2k3  \ 

- s 2  -  k2 
Sa/s2  +  \k  2 

\ 

( 

{ 

\  s^/s2  +  \k\2  ) 

V  0  ) 

J 

Straightforward  calculations  show  that  for  this  subspace  the  solution  components 
obey  the  relationships 

(2.8) 


, _  7.2  1 

\/s2  +  \k\2  +  s  +  -3- 
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s  ,  k2k3- 

-C/2 - ~ — c/3  H - ~ — 2  — 


_ _  b2~\ 
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Bt  =  0, 


(2.9) 


-^4  + 


\/s2  +  \k\2  +  s  + 


E3  + 


\fs2  +  \k\2  +  s  + 


b2  -  =  o. 

s 


The  earlier  review  article  [35]  discusses  the  origin  of  these  relationships  in  terms  of 
the  left  eigenvectors  of  the  matrix  appearing  in  (2.5). 

To  produce  physical-space  time-domain  radiation  conditions  from  (2.8)-(2.9), 
we  must  carry  out  the  necessary  inverse  transformations.  We  first  introduce  the 
kernel  [35,  3] 

(2.10)  K(t)  =  t~1J1(t),  K(s)  =  Vs2  +  l-s, 
where  of  course  K(s)  decays  in  s.  In  (2.8)-(2.9)  we  set 

(2.11)  y/§2TW  +  s=\k\K(\k\-1s)+2s, 

rearrange  terms,  and  make  substitutions  with  (2.4),  in  order  to  reach  a  set  of 
equations  on  which  the  inverse  transformations  are  easily  carried  out.  We  then  find 

(2.12)  H®(j*-B3)  +  W*-B3)  +  ! 

(2.U)  lJt(E3  +  B2)  +  K{E3  +  B2)+aJ±  +  aJh.  =  o, 

where  the  nonlocal  operation  (IZw)  (0,  y,  z,  t)  is  defined  through 

rt  j  ('clfclr) 

(2.14)  J7(V,w)(0,k2,k3,t)  =  /  - — - [c\k\2w{0,k2,k3,t-T)]dT. 

Jo  c\k\t 

Here  T  denotes  Fourier  transform  in  the  tangential  variables  (y,  z ),  and  succinctly 

(2.15)  Uw  =  T~x  (c\k\2K(c\k\t)  *  {Fw)) . 

We  note  that  analogous  expressions  can  be  derived  in  waveguides,  which  is  the 
most  practical  application  of  the  planar  boundary  formulas.  We  consider  only  the 
simplest  possible  case.  In  particular,  suppose  that  for  x  >  0  the  waveguide  has 
constant  rectangular  cross-section,  0  =  [0, x  [0 ,LZ].  Suppose  further  that  the 
walls  are  perfectly  conducting.  Then  we  may  replace  the  Fourier  transfroms  in  the 
expressions  above  by  Fourier  series.  (Note  that  for  more  complicated  cross-sections 
the  relevant  eigenfunction  expansions  couple  Cartesian  components  in  a  nontrivial 
way.)  Noting  the  boundary  conditions  and  divergence  constraint  we  conclude  that 
E2  and  B3  should  be  expanded  in  terms  of  cos  (^- 77-^  •  sin  j  =  CSk2,k3  while 
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E3  and  B2  are  expressed  in  terms  of  sin  •  cos  =  SCk2,k3 ■  (See,  e.g, 

[61,  Ch.  9].)  We  must  now  only  replace  T  by  the  series  transform  Fes  in  equation 
(2.12)  and  Fsc  in  equation  (2.13),  noting  that  now  |fc|2  =  7r2  . 

Refs.  [3,  50]  present  efficient  strategies  for  numerically  implementing  the  con¬ 
volution  (2.15),  and  these  methods  will  be  discussed  below.  Such  approximations 
obviate  the  need  to  carry  out  the  exact  inverse  Fourier  transform  in  (2.15).  Nev¬ 
ertheless,  as  an  interesting  exercise  we  here  perform  this  inverse  transformation  in 
order  to  achieve  the  exact  physical-space  time-domain  boundary  conditions.  First, 
with  y  =  (y,  z),  k  =  (k2,k3),  u  =  (u,v),  and 

(2.16)  F(k  2,k3)  =  ^-[  e~lykf(y,  z)dy,  f(y,z)  =  -^-[  e*y'kF(/c2,  k3)dk, 

J R2  2tt  Jr 2 


recall  that  by  the  Fourier  convolution  theorem 

(2.17)  F~1[F(k2,k3)G(k2,k3)](y,z)  =  T  j  f(y-u,z-v)g(u,v)du. 

47T  J R2 


In  order  to  apply  the  convolution  theorem  (2.17),  we  first  assemble  several  results 
from  Watson’s  monograph  [74]  in  order  to  compute  the  inverse  Fourier  transform 


(2.18) 


J_  [  iykMc\k\t)  =  H{t-  p/c) 
2n  2  c\k\t  c2t2 


where  p  =  \Jy2  +  z2  and  i/(£)  is  the  Heaviside  step  function  such  that  H{£)  =  0 
for  £  <  0,  17(0)  =  i?(£)  =  1  for  £  >  0.  With  these  results,  we  find  that 

(2.19)  (1lw)(0,y,z,t)  =  -2-  f  1  f  Auu>(0,  u,v,t )dudr , 

2tt  J o  c(t  -  tY  J|y-u|<c(t-T) 

with  Au  denoting  the  Laplacian  in  the  u-v  plane.  The  divergence  theorem  then 
gives 


(lZw)(0,  y,  z,  t)  = 

1  F  l  d  f2n 

(2.20)  —  —  /  — - -  —  /  u>(0,  y  +  c(t  —  t)  cos  9,  z  +  c(t  —  r)  sin  9,  r)d0dT. 

2tt  J0  c(t  -  t)  at  J 0 

Notice  that  the  integration  is  over  the  intersection  of  the  artificial  boundary  x  =  0 
and  the  past  liglitcone  belonging  to  the  spacetime  point  (0,y,z,t). 


2.2.  Spherical  boundary. 


2.2.1.  Vector  spherical  harmonics.  We  consider  both  pure-spin  and  pure-orbital 
vector  spherical  harmonics  (see  [69]  for  the  origins  of  this  terminology).  As  given 
in  [35,  53,  56,  55],  the  unnornralized  pure-spin  harmonics  are  the  set 


(2.21) 


Y  Pm  — Vprrl  er , 

lTf  JJYtrn  ,  1  dYlm 

^  im  —  p.Q  e0  H  : — X  o  j.  e 
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sin  0  d(p 
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where  er,  eg,  and  are  the  standard  unit  basis  vectors  in  the  spherical  polar 
coordinate  system  and 


(2.22) 


Ytm(9,4>)  = 


4t  (£  +  to)  ! 


eim<t>p™  (cog  0), 


is  one  of  the  standard  spherical  harmonics,  orthogonal  on  the  unit  sphere.  (Here 
P ™  is  the  associated  Legendre  function  defined  in  [1],  Cli.  8.)  When  paired  with 
r-dependent  expansion  coefficients,  these  harmonics  are  easily  seen  to  form  a  closed 
system  under  the  standard  vector  operations  (div,  grad,  curl)  involving  the  gradient 
operator  V  [55].  In  terms  of  the  pure-spin  harmonics,  we  also  define  a  set  of 
normalized  pure-orbital  vector  harmonics, 


=Yt1 


(2.23)  X,m  =Y|m  =  —I 
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^  VP/; 
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V  2£+l 

ww+ 1)  J 
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J  e 
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The  Y[m  are  Thorne’s  Y^  >em  [69],  and  the  V(m,  X^TO,  W (m  notation  is  due  to 
Hill  [47].  When  paired  with  r-dependent  expansion  coefficients,  the  pure-orbital 
harmonics  also  form  a  closed  system  under  the  standard  vector  operations  involving 
the  gradient  operator  V,  and  a  compendium  of  formulas  is  given  by  Hill  [47].  Apart 
from  the  factors  of  1/ \J  £{£  +  1),  which  serve  to  normalize  the  pure-spin  harmonics 
(2.21),  the  matrix  associated  with  the  transformation  (2.23)  is  unitary. 

The  orbital  harmonics  (2.23)  also  arise  in  the  quantum  theory  of  angular  mo¬ 
menta  [11].  Indeed,  in  terms  of  the  complexified  basis  to  =  ez,  t-j-i  =  ^(e^,  A 
iey  )/v/ 2,  they  are 


f  l 

(2.24)  Y(m  =  J2  Y 

m'=— i'  m"=— 1 


where  the  (£' \m' m" \P l£m)  are  Clebsch-Gordan  coefficients  [1,  11].  For  fixed  t 
the  Y|m  transform  under  an  order-f  representation  of  the  rotation  group  SO( 3). 
The  expansion  above  expresses  this  representation  as  a  coupling  between  the  scalar 
harmonics  Y^iTO  (an  order-!?'  representation,  where  £'  =  £ —  1,£,£  +  1)  and  the  basis 
tm  (an  order-1  representation).  Starting  from  (2.23)  and  using  identities  for  the 
scalar  harmonics  (collected,  for  example,  in  [48]),  one  may  also  directly  calculate 
the  following  explicit  expressions: 


t~\  (£  +  m)(£  +  m- 1) 


vt-i  _ 

1  Cm  — 


2£{2£  —  1) 
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(2.25) 
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The  factors  involving  square  roots  are  the  nonzero  Clebsch-Gordan  coefficients  in 
the  expansion  (2.24).  Unlike  the  pure-spin  harmonics  (2.21),  the  Y|m  are  eigen¬ 
functions  of  the  spherical  Laplacian  with  eigenvalue  +  1). 


2.2.2.  Exact  radiation  boundary  conditions.  Using  the  pure-spin  harmonics,  Ref.  [35] 
has  formulated  exact  radiation  boundary  conditions  for  the  electromagnetic  held  in 
the  presence  of  a  spherical  boundary.  Here  we  will  provide  an  equivalent  description 
in  terms  of  the  pure-orbital  harmonics.  Nevertheless,  since  they  are  tailored  to  the 
transverse  character  of  the  radiation  held,  we  at  first  work  with  the  pure-spin  har¬ 
monics,  only  converting  to  the  pure-orbital  harmonics  once  our  calculations  have 
been  completed. 

To  start,  we  perform  a  Laplace  transform  on  the  system  (1.1)— (1.2),  subsequently 
expanding  the  transformed  variables  E  and  B  in  pure-spin  harmonics.  The  E 
expansion,  for  example,  is 

OO  £ 

(2.26)  =  +  E™ *Cm  . 

£—1  m=—£ 


This  process  leads  to  a  differential-algebraic  system  of  equations,  where  the  alge¬ 
braic  sector  is 


(2.27) 


~sErtm  =  - 


'+1) 
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(2) 
£m  ’ 


SBrtm  = 


l± 

r  ^ Cm ' 


These  equations  may  be  used  to  eliminate  the  radial-harmonic  coefficients,  E\m 
and  Bpm,  and  such  an  elimination  yields  the  following  first-order  system  for  the 
remaining  coefficients: 

(2.28) 
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The  solutions 
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remain  finite  as  r  — >  oo  for  Re(s)  >  0.  Here  kg(z)  is  a  modified  spherical  Bessel 
function  expressed  as  kt(z)  =  \J n / (2z)K i2(z)  in  terms  of  the  standard  modi¬ 
fied  cylindrical  Bessel  function  Kv(z)  known  as  the  MacDonald  function  [74].  As 
defined, 

(2.30)  ke(z)  ~  7r(2z)~1e~z  as  z  — >  oo, 

but  some  authors  fix  the  definition  of  kt{z)  so  that  kg(z)  ~  z~1e~z.  In  any  case,  the 
choice  of  the  overall  of  overall  constant  does  not  affect  our  argument.  From  now  on 
we  assume  that  the  multipole  components  have  the  stated  form,  so,  for  example, 
=  ~aim{s)kf,{sr).  Further  calculations  then  show  that 

(2-31)  +  B%)  +  -rm~sr)E%  =  0,  3(E%  -  B<%)  -  ^Me(Sr)B^  =  0. 

Here  the  frequency-domain  kernel  is 

1  _^+i/2(z) 

2  K(+1/2{z) 

In  fact  M((z)  =  0(1/ z),  so  that  overall  r_1Mp(sr)  is  0(l/r2)  in  the  radial  coordi¬ 
nate. 

By  the  Laplace  convolution  theorem,  the  corresponding  time— domain  boundary 
conditions  are  then  (with  r  =  R  taken  as  the  boundary  location) 

(2-33)  +  B<»)  +  I  J*  Me(cr/R)E2(R,t  -  r)dr  =  0, 

(2'34)  -  Btm)  ~  Ji  l  Me(cT/R)B<V(R,t-  r)dr  =  0, 

where  the  time-domain  kernel  is 


(2.32) 


Me(z)  =  - 


1  +  z  +  z 


KK) 

kAz)  I 


t 

(2.35)  M((ct/R)  =  -  ^2(zetk/R)  exp (zi^ct/R) 

k= i 


in  terms  of  the  t  roots  zetk  of  Ki+1/2(z)  which  all  he  in  the  left-half  plane  2-plane. 

As  will  be  discussed  in  subsequent  sections,  the  fact  that  the  kernels  are  expo¬ 
nential  functions  of  t  implies  that  the  temporal  convolution  can  be  localized.  That 
is,  one  can  avoid  the  storage  of  the  time  history  evident  in  (2.33)-(2.34).  This  fact 
was  independently  exploited  by  Sofronov  [63,  64]  and  Grote  and  Keller  [28,  29,  30] 
to  derive  and  implement  accurate  temporally  local  conditions.  (See  also  [31]  for 
applications  to  multiple  scattering.)  In  fact,  in  [38]  it  is  shown  that  purely  lo¬ 
cal  conditions  can  be  developed  which  are  exact  for  solutions  described  by  finitely 
many  spherical  harmonics.  In  particular,  setting  wq  =  2 u,  rop+i  =  0  and  solving 
the  following  coupled  sequence  of  equations  on  the  sphere: 

1  dwj  j  1  .. . 

-~Qf  +  -rw3  =  ^2  (As2  +3V  -  !))  w3- 1  +  wj+ 1.  j  =  !>  ■  •  • >  K 


(2.36) 
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it  is  shown  that: 

p  i 

(2.37)  W!  +  -  EE  Yim  [Ml  *  Uem)  =  0, 

£—0  m=—£ 


if  we  assume  uem  =  0  for  l  >  P.  This  local  exact  condition  has  been  successfully 
implemented  by  Grote  [27]. 

These  formulations  become  expensive  as  the  harmonic  index,  £,  increases  as 
(2.35)  requires  i  exponentials.  However,  as  described  in  [2],  for  large  (,  one  may  use 
compressed  kernels  and  fewer  exponential  terms  while  retaining  high  accuracy.  We 
will  describe  this  development  below. 

Introducing  B  =  erxB  =  -B^eg+Bge^,  so  that  B^  =  and  -B^  =  B^, 
we  then  write  (2.33)-(2.34)  as  follows: 

(2-38)  +  E%)  =  0 

(2-39)  +  ^Mt  *  =  °- 

Finally,  summing  these  equations  on  the  harmonics,  we  obtain 


(2.40) 


-  —  (  Ft 

cdv 


5)  +  ^E  E 

£=1  rn=—£ 


$ 


£m 


where  the  superscript  T  denotes  “transverse”,  that  is  ET  =  E  —  E(E  ■  er).  Note 
that  Et  +  B  has  components  (0,  Eg  —  B$,  Eg  +  B^). 

We  now  rewrite  (2.40)  in  terms  of  the  pure-orbital  harmonics,  thereby  achieving 
an  expression  which  may  be  implemented  using  the  scalar-harmonic  transform. 
Since  B  =  BT ,  it  follows  from  (2.23)  that 


(2.41)  C_1)  =  +  VBtZ  C+1)  = 

Here  the  expansion  coefficients  are  with  respect  the  pure-orbital  harmonics, 


(2.42) 


=  /  B  ■  Ye;mdS, 
Js2 


where  the  overline  on  Y  indicates  complex  conjugation  and  dS  is  the  area  measure 
on  the  unit  sphere  S2.  In  terms  of  the  pure-orbital  harmonics,  the  boundary 
condition  (2.40)  becomes 


Id,  T 


B) 

E  E  [ 

£=1  m——£ 


•y£— 1 
1  £m 


(2.43) 


Taking  advantage  of  the  list  (2.25),  we  may  use  the  scalar-harmonic  transform 
to  compute  the  coefficients  B^J  and  E^J  .  This  fact  is  significant  for  practical 
implementations  as  we  can  use  well-developed  software  for  computing  and  inverting 
the  transform.  (See,  for  example,  the  routines  distributed  with  the  NCAR  spectral 
transform  shallow  water  model  [57].) 
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2.3.  Cylindrical  boundary.  Our  last  example  is  a  cylindrical  boundary  of  infi¬ 
nite  extent,  defined  in  terms  of  the  standard  cylindrical  coordinate  system  (r,  9 ,  z) 
as  r  =  R.  We  expand  the  vector  field  E  in  terms  of  the  standard  orthonormal 
cylindrical  basis  {er,e@,ez}  as  E  =  Erer  +  Egeg  +  Ezez,  and  similarly  for  B. 
Next,  we  write  down  the  Maxwell  system  (1.1) — (1.2)  component-by-component, 
subsequently  performing  on  each  field  component  a  Laplace  transform  in  t  (with 
dual  variable  s  and  denoted  by  a  hat),  continuous  Fourier  transform  in  z  (with 
dual  variable  k  and  denoted  by  a  bar),  and  a  Fourier  series  expansion  in  9  (with 
dual  index  n  and  denoted  by  a  superscript  n).  Like  before,  this  process  yields  both 
algebraic  equations, 


(2.44)  sE?  -  B nz  -  ikB^j  =  0,  SB"  +  (^jEnz  -  ikE^j  =  0, 

as  well  as  a  system  of  ordinary  differential  equations, 

(2.45) 


/  rEg 

\ 

(  0 

0 

_ kn 

sr  „ 

Sr+£  \ 

1  sr 

( rk 

\ 

k 

0 

0 

_i  _  j£_ 

kn 

k 

z 

rk 

+ 

kn 

sr 

—sr  —  k- 

sr 

r  sr 

0 

sr 

0 

z 

rk 

\  k 

) 

\  l  +  ±- 
\  r  sr 

kn 

sr 

0 

o  ) 

\  k 

/ 

The  solutions  which  remain  finite  as  r  — >  oo  have  the  form 


(2.46) 


/  rEg  \ 

k 

rBng 

\  k  ) 


/ 

=  an(s,  k) 

V 


%Kn{ra) 
Kn  (ra ) 
f  Kira) 
0 


)  ( 

+  bn(s,  k) 


J  V 


-fK(ra)  \ 

0 

%Kn(r<r) 
Kn(ra)  J 


where  cr  =  \/s2  +  k2,  with  the  branch  chosen  as  in  plane  boundary  example.  Notice 
that  cr  —  S  =  kK (fc-1S),  where  the  kernel  K(s),  first  appearing  in  (2.10),  should  not 
be  confused  with  the  MacDonald  function  Kn(ra).  For  the  chosen  subspace  (2.46) 
of  solutions  the  components  satisfy  the  relationships 


(2.47)  s(k  +  k)  +  Y~+  ikk  +  kkik-^k  +  ~Cn  ( rVs 2  +  k2)Enz 

Brl 


(2.48)  S(.B"  -  Eg)  +  +ikB ™  +  kk{k~1s)Bnz  +  -Cn{r^/~s2  +  k2)B"z 

To  reach  these  equations  we  have  made  use  of  (2.44)  and  introduced 


(2.49) 


Cn(z)  =  - 


1 

2+2’ 


K(z) 


Kn(z)  \ 


=  0, 
=  0. 


Since  Cn(z )  ~  z~1  as  z  —>  oo,  the  kernel  Cn(rV s 2  +  k2)  is  the  Laplace  trans¬ 
form  of  a  function  Gn(ct,r,k).  We  may  now  easily  perform  the  requisite  inverse 
transformations  on  (2.47)-(2.48).  Evaluated  at  the  boundary  r  =  R,  the  results 
are 

-^-(Ez  +  Bg)  +  k  +  ^L+7ZzEz  +  QgzEz  =  0, 
cot.  2  R  oz 

l^-(Bz-Eg)+^-  +  kk  +  KzBz  +  QgzBz=  0. 
cot  2R  oz 


(2.50) 

(2.51) 


10 


THOMAS  HAGSTROM  AND  STEPHEN  LAU 


Here,  with  Tz  indicating  Fourier  transform  in  2  and  w  =  w(R ,  91  2,  t),  we  define 

(2.52)  Kzw  =  T~x  (ck2K(ckt)  *  {Fzw)) , 

which  is  similar  to  formula  (2.15).  Moreover,  with  Tgz  representing  double  Fourier 
transformation  in  9  (series)  and  2  (continuous),  the  remaining  nonlocal  operation 
is 

(2.53)  Qgzw  =  TgX  ( Gn(ct ,  R,  k)  *  Rgzw) . 

Although  we  will  not  explicitly  perform  the  inverse  Laplace  and  Fourier  transforms 
in  order  to  define  the  exact  time-domain  physical-space  boundary  conditions,  we 
note  that  efficient  strategies  exist  for  numerically  implementing  these  nonlocal  op¬ 
erations  [3]. 

2.4.  Fast  time-local  evaluation  of  the  kernels  I:  global  exponential  ap¬ 
proximations.  Clearly,  a  primary  bottleneck  in  the  direct  evaluation  of  the  exact 
boundary  conditions  derived  above  is  the  evaluation  of  integral  operators.  To  give 
a  rough  count  of  the  complexity,  we  suppose  lengths  are  scaled  by  the  dimensions  of 
the  computational  domain,  and  time  is  scaled  so  that  c  =  1.  Then  the  integration 
time,  T,  measures  the  number  of  times  a  wave  could  traverse  the  domain.  If  a 
characteristic  wavelength  in  these  units  is  given  by  A  the  complexity  of  a  standard 
solver  would  be: 

(2.54)  Work  oc  A "4T,  Storage  a  A“3. 

A  reasonable  goal  is  that  the  cost  of  boundary  treatment  is  no  worse  than  these. 
Specializing  to  the  spherical  boundary,  the  cost  per  time  step  associated  with  the 
integrals  is  the  sum  of  the  cost  of  the  spherical  harmonic  transform  and  the  cost 
of  the  temporal  convolution.  Noting  that  the  number  of  harmonics  required  will 
scale  like  A-2,  the  use  of  standard  software  such  as  [57],  which  combines  FFTs  in 
the  azimuthal  coordinate  with  direct  transforms  in  6,  leads  to  a  total  cost  of: 

(2.55)  WorksHT  ex  A"4T, 

which  is  comparable  to  (and  in  practice  less  than)  the  work  required  by  the  volume 
solver.  For  A  <C  1  one  could  instead  use  one  of  the  recently  developed  fast  spherical 
harmonic  transforms  (e.g.  [52,  43,  65]).  These  reduce  the  complexity  to: 

(2.56)  WorksHT  cx  A“3  In  (A_1)T. 

The  work  associated  with  the  temporal  convolution  can  also  be  kept  manage¬ 
able  through  the  use  of  fast  alogrithms.  Precisely,  Hairer  et  al  [41]  show  that  the 
well-known  FFT-based  algorithms  for  computing  convolutions  can  be  adapted  to 
evolutional  convolution  integrals  of  our  form.  Using  their  algorithm  we  have: 

(2.57)  Workconv  oc  A“3  In2  (A ~X)T  In2  T. 

However,  the  algorithm  in  [41]  requires  full  storage  of  the  boundary  history: 

(2.58)  Storageconv  oc  A _3T. 

This  is  the  dominant  storage  cost  for  T  large,  and  is  prohibitive  for  large  applica¬ 
tions. 
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The  key  idea  in  developing  fast,  low-storage  implementations  of  the  exact  bound¬ 
ary  conditions  is  the  observation  that  convolution  with  exponential  functions  re¬ 
quires  no  history  storage.  Precisely,  if: 

(2.59)  <j>{t)  =  f  ae/3(t~r)u(r)dT, 

Jo 

then  <j>  satisfies  the  differential  equation: 

(2.60)  =  p(f>  +  au. 

at 

We  can  thus  expect  to  compute  <f>  using  0(X~1T)  work  and  0(1)  storage,  vastly 
improving  on  the  general  algorithm  of  [41].  For  the  case  of  a  spherical  boundary, 
we  have  already  observed  that  the  kernels  are  exactly  equal  to  sums  of  exponentials 
(2.35).  This  leads  to  an  algorithm  which  is  mathematically  equivalent  to  the  one 
proposed  by  Sofronov  [63,  64]  and  Grote-Keller  [28,  29,  30]: 

(2-61)  f  Me(c(t  -  t)/R)  (  ^  TJ  dr  =  £ 

J 0  V  Btm(R’T)  J 


(2.62) 


dZe,k 

dt 


Z£,k,C  rj  Z£,kC 

-£*etk  — 


R 


R 


p(2) 

J  im 


e)2{r,t) 


Z/,fc(0)  =  0. 


Now  the  number  of  auxiliary  functions  Z e,k{t)  which  must  be  computed  is  propor¬ 
tional  to  A-3.  Thus  the  cost  of  the  local  convolution  algorithm  is: 


(2.63) 


Workc 


oc  A  4T,  Storageconv  oc  A 


which  is  comparable  to  the  work  and  storage  required  by  the  volume  solver. 

Further  improvements  in  efficiency  and  generalizations  to  the  planar  and  cylindri¬ 
cal  boundaries  follow  from  the  uniform  approximation  of  the  temporal  convolution 
kernels  by  a  smaller  number  of  exponentials.  The  analysis  and  practical  construc¬ 
tion  of  these  approximations  is  carried  out  in  [2,  3].  We  will  review  their  analysis 
briefly.  However,  from  the  user’s  perspective,  one  only  needs  a  table  of  exponents 
and  amplitudes.  These  can  be  obtained  at  [34]. 

Let  Cp(t)  denote  any  of  the  convolution  kernels  derived  above  (see  (2. 15), (2. 43), (2. 52)) 
with  p  representing  the  spatial  harmonic  index.  The  approximation  problem  is  to 
find  (ctj'p,  (3jjP),  j  =  1,. . .  ,NP,  such  that: 

Np 

(2.64)  Rp(t)  =  Y/aj,pe/3^t, 

3=1 

satisfies  for  some  (small)  tolerance  e  and  (large)  time  T: 

(2-65)  || (RP  -  Cp)  *  /||l2(o,t)  <  ell/IU2(o,T). 

By  Parseval’s  relation  we  can  translate  this  to  an  equivalent  statement  on  the 
rational  approximation  of  the  Laplace  transform  of  the  kernel: 

max  | Rv ( s )  —  CVAs)!  < 

3?(s)=t-1  F  F  e 


(2.66) 
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where  we  note  that: 

Np 

(2.67)  ^)  =  E7^- 

j-i  *  PAp 

The  fundamental  theoretical  result  of  [2,  3]  is  that  all  the  kernels  admit  expo¬ 
nential  approximations  with: 

(2.68)  Np  <  C In  -  •  In  ( pT ). 

The  proof  is  based  on  representations  of  Cp  as  sums/integrals  of  poles  over  ap¬ 
propriate  contours  in  the  left  half  s-plane  combined  with  a  general  approximation 
result  for  functions  so  represented.  Using  these  exponential  approximations  the 
convolution  algorithm  now  costs: 

(2.69)  Workconv  a  A~3Tln  —  ■  In  >  Storageconv  a  A-2  In  -  •  In  (^'j  ■ 

Now  the  work  and  storage  is  negligible  in  comparison  with  the  requirements  of  the 
volume  solver.  Only  the  cost  of  the  spatial  transforms  is  formally  comparable. 

The  practical  numerical  construction  of  the  poles  and  amplitudes  yields  remark¬ 
ably  efficient  exponential  approximations.  In  Table  1  we  list  results  for  a  tolerance 
of  e  =  10~6.  (The  poles  and  amplitudes  themselves  are  available  at  [34].)  Note  that 
we  exploit  the  homogeneity  of  the  plane  kernel  to  compute  fc-independent  poles: 

(2.70)  ccj.fe  =  \k\2aj,  (3j}k  =  \k\0j. 

Also  note  that  the  zero  mode  for  the  cylinder  kernel  is  particularly  difficult  to 
approximate,  but  the  highers  modes  are  essentially  the  same  as  in  the  spherical 
case. 


Cylinder 

Sphere 

Plane 

n 

Nn 

i 

Nt 

\k\T  <  104 

0 

26 

31 

1 

9 

2 

6 

3-6 

5 

0-5 

l 

7-8 

6 

6-8 

6 

9-12 

7 

9-12 

7 

13  -  19 

8 

13-19 

8 

20-31 

9 

20-31 

9 

32-51 

10 

32-51 

10 

52-86 

11 

52-86 

11 

87-147 

12 

87  -  147 

12 

148  -  227 

13 

148  -  228 

13 

228  -  401 

14 

229  -  402 

14 

402  -  728 

15 

403  -  728 

15 

729  -  1024 

16 

729  -  1024 

16 

Table  1.  Poles  required  for  exponential  approximations  to  the 
nonreflecting  boundary  kernels  with  a  tolerance  of  e  =  10~6. 
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2.5.  Fast  time-local  evaluation  of  the  kernels  II:  piecewise  exponential 
approximations.  More  recently,  alternative  fast,  low-storage  algorithms  for  eval¬ 
uating  evolutional  convolutions  have  been  proposed  by  Schadle  and  coworkers 
[50,  49,  60].  These  are  based  on  piecewise  rather  than  global  exponential  approxi¬ 
mations  to  the  kernel.  The  advantages  of  using  piecewise  rather  than  global  approx¬ 
imations  are  that  they  are  easier  to  construct  and  more  generally  applicable.  The 
disadvantages  are  that  the  very  simple  convolution  algorithm  embodied  in  (2.60) 
must  be  replaced  by  an  algorithm  which  is  more  complicated,  and  slightly  more 
memory  and  work  are  required. 

Concerning  generality,  it  is  shown  in  [60]  that  the  only  essential  requirement  for 
the  applicability  of  the  approach  is  that  the  Laplace  transform  of  the  convolution 
kernel,  C(t),  be  sectorial.  That  is,  for  some  complex  number  sq  and  angle  </>  < 

C(s )  is  analytic  for  |arg(s  —  So)|  <  tt  —  </>.  Moreover,  for  positive  constants  M  and 
v. 

(2.71)  \C(s)\  <  M\s\~v. 

These  requirements  are  satisfied  by  the  kernels  Cp(t)  discussed  above,  but,  of  course, 
by  many  other  kernels  appearing  in  diverse  applications.  The  exponential  approx¬ 
imations  themselves  are  defined  on  intervals  of  rapidly  increasing  length.  Fixing  a 
base,  B  >  1,  and  a  time  step  At  the  kernel  is  approximated  by  a  fixed  number  of 
exponentials,  P  oc  In  on  each  subinterval: 

p 

(2.72)  C{t)  w  t  €  [5fc_1  At,  (2 Bk  -  l)Ai]  =  Ik. 

i=i 

(Notice  the  overlap.)  The  poles  and  amplitudes  are  directly  computed  from  C(s) 
by  applying  a  P- point  quadrature  rule  to  the  Laplace  inversion  integral  along  a 
specially  chosen  contour;  the  so-called  Talbot  contour  [66].  For  fixed  B  one  can 
take  P  oc  In  -  for  an  error  tolerance  e.  Thus  the  least  squares  procedure  of  [2,  3] 
is  avoided,  and  the  approximations  can  be  computed  on  the  fly;  all  the  user  needs 
to  know  is  the  singularity  structure  of  C  which  determines  some  parameters  in  the 
contours. 

The  approximate  convolution  may  be  derived  as  follows.  If  the  final  simulation 
time  is  T  it  is  clear  that  the  piecewise  approximations  to  C  will  eventually  be 
needed  on  intervals  Ik,  k  =  1  where  L  is  the  smallest  integer  such  that 

(2 Bl  —  l)Af  >  T.  Clearly  L  oc  ln(A_1T).  We  solve  the  differential  equations 
associated  with  (2.72)  over  intervals  ((£  —  l)BkAt,£BkAt),  £  =  1,...,£f,  where 
(£f  +  l)BkAt  >  T.  Precisely,  suppose 

(2J3)  =  pf- +  af yjikit  ((£  -  l)Bk  At)  =  0, 

and  set: 

(2.74)  zjtktttP  =  yjtkit(((e-l)Bk+pBk-1)to),  p  =  l,...,B. 

Note  that  the  work  involved  in  computing  these  numbers  is  propoprtional  to  PLX~1T  oc 
Ini  .  In (A-1T)A-1T.  Now  consider  the  approximate  evaluation  of 

(2.75)  f  C(t  —  t)u(t)cIt. 

Jo 


14 


THOMAS  HAGSTROM  AND  STEPHEN  LAU 


Obviously  the  local  part  of  the  integral  can  be  approximated  using  local  data  and 
so  we  concentrate  on  the  integral  up  to  t  —  At.  We  partition  it  into  subintervals 
over  which  the  various  kernel  approximations  are  valid.  Precisely,  we  determine 
t  —  At  =  To  >  Ti  >  ...  >  tm  =  0  such  that  Tk  =  {Ik  —  l)Bk At  for  some  integer  Ik 
and  ( t  —  Tk- 1 ,  t  —  Tk)  G  Ik-  (The  construction  of  the  partition  involves  an  expansion 
of  the  time  step  index  in  base  B;  see  [49].)  We  then  have: 

pt-At 

/  C(t  —  t)u{t)cLt  = 

J  o 

(2.76)  « 

where  {Ik  —  1  )B +pk  =  Ik- 1  —  1.  Applying  this  algorithm  to  the  nonlocal  boundary 
conditions  described  above  and  being  careful  about  which  values  Zj^k,i,p  can  be 
discarded  we  find  that  the  cost  is: 

(2.77)  Workconv  oc  A~3Tln  -  •  In  (A_1T),  Storageconv  oc  A-2  In  -  •  In  (A-1T), 

which  is  again  negligible  in  comparison  with  the  volume  solver. 

Of  course,  it  can  be  argued  that  these  developments  do  not  improve  on  the 
existing  global  approximations  which  have  already  been  constructed  and  tabulated, 
and  which  we  have  seen  to  be  quite  efficient  for  the  kernels  arising  in  applications 
to  Maxwell’s  equations.  However,  they  can  be  used  to  evaluate  exact  conditions  for 
spatial  discretizations,  which  may  prove  more  accurate  for  problems  with  unresolved 
waves  of  nonnegligible  amplitude.  This  is  carried  out  in  detail  in  [49].  Furthermore, 
they  may  be  useful  for  generalizations  to  more  complex  media. 


E  /  C{t~T)u{T)dT 
fc=  1  "'Tfc 


EE' 

fc= i  i= i 


ePjk)(t-Tk- i) 

e  ^j,k,tk,Pk 


3.  Exact  formulations  on  general  boundaries 

Despite  the  existence  of  fast,  low-storage  evaluation  algorithms,  the  fact  that  the 
boundary  conditions  considered  so  far  require  the  use  of  a  restricted  set  of  artificial 
boundaries  does  lead  to  nonnegligible  costs  in  some  cases.  In  particular,  the  need 
to  embed  a  high-aspect-ratio  scatterer  in  a  spherical  computational  domain  may 
drastically  increase  the  computational  requirements.  Therefore  one  would  like  to 
construct  accurate  radiation  conditions  on  more  general  boundaries,  and  we  will 
discuss  such  procedures  for  the  remainder  of  the  article. 


3.1.  Kirchhoff  representations.  Rather  early  on  Ting  and  Miksis  [70]  proposed 
a  scheme  for  implementing  exact  boundary  conditions  for  the  time-dependent  scalar 
wave  equation.  They  considered  a  scenario  as  in  Figure  1,  with  the  computational 
domain  extended  spatially  and  taken  to  lie  within  another  artificial  boundary  T'. 
Provided  that  the  initial  data,  w(-,0)  and  itt(-,  0),  lies  within  T  (and  further  that 
any  inhomogeneities  are  confined  both  spatially  and  temporally  to  the  region  of  H 
within  r),  then  field  values  on  T'  have  a  retarded  Kirchhoff  representation,  [9] 


(3.78) 
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where  dSx  is  the  area  measure  and  n  is  the  inward-pointing  normal  on  F  (that  is 
outward-  pointing  with  respect  to  the  tail  H).  Givoli  and  Kohen  [25]  numerically 
implemented  this  scheme  for  both  the  scalar  wave  equation  in  three  space  dimen¬ 
sions  and  for  the  equations  of  elasticity.  He  and  Weston  [42]  developed  a  fully 
vector  version  of  the  scheme  as  applied  to  the  Maxwell  equations. 


Figure  1.  Domains  for  an  exterior  problem.  H  is  the  tail, 
Q  is  the  computational  domain,  and  both  T  and  T'  are  artificial 
boundaries.  The  irregular  objects  within  O  represent  scatterers. 


Teng  [68]  demonstrated  that  one  can  do  away  with  the  need  for  two  artificial 
boundaries,  in  effect  considering  the  limit  when  T  and  T'  are  the  same  surface.  For 
a  generic  point  x'  G  T,  he  finds  that 


i(x',t)  1  - 


Q(x') 

47T 


1 

47T 


i(x,  t  —  r/c) 


dn 


(3.79) 


1  du ,  .  1  dr  du ,  .  . 

“i-  x.^r/c - o-  x-^r/c 

r  on  rc  on  at 


dS x 


x'  g  r, 


where  B(x')  is  the  exterior  solid  angle  as  measure  in  the  tangent  space  at  x'. 
Provided  that  the  artificial  boundary  is  smooth  at  x',  his  formula  reduces  to 


T  L 


/  ,  .  d  ( 1 

„(*,(-,•/<%  i- 


(3.80) 


1  du  .  .  1  dr  du  .  . 

- -5-(x>*-  r  c) - r/c) 

r  dn  rc  dn  dt 


dSx 


x'  g  r. 


This  boundary  condition  is  clearly  nonlocal  in  space  and  time;  however,  its  history 
dependence  is  restricted  in  the  following  sense.  The  integral  involves  the  retarded 
time  t  =  t  —  r/c  which  may  be  confined  to  the  interval  t  —  rmax/c  <  r  <  t,  where 
r max  is  the  maximum  Euclidean  distance  between  any  two  points  on  T. 


3.2.  Origin  of  simple  of  history  terms.  Teng’s  derivation  of  (3.79)  relies  on  the 
theory  of  distributions.  Although  we  shall  not  repeat  his  argument,  let  us  show  how 
such  history-dependent  terms  arise  in  a  restricted  setting,  that  is  the  homogeneous 
scalar  wave  equation  for  c  =  1  and  a  simple  class  of  infinite-extent  boundaries, 
such  as  an  infinite  plane,  two  semi-infinite  planes  which  meet  at  an  edge,  or  three 
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quarter-infinite  planes  which  meet  at  a  corner.  We  will  derive  a  formula  for  w(x',  T) 
at  time  T,  with  the  spatial  point  x'  =  (0,0,0)  taken  as  the  coordinate  origin  and 
assumed  common  to  all  the  planes  which  make  up  the  boundary  F.  Let  B(r,t) 
represent  the  radius-r  sphere  at  time  t  and  centered  at  the  origin.  With  this 
notation  we  write  B[T  —  t,t)  for  the  intersection  of  time  level  t  <  T  and  the 
past  lightcone  of  the  spacetime  point  (0,  0,0,  T).  The  artificial  boundary  T  divides 
B(T  —  t,t)  into  two  components,  each  one  a  spherical  polygon,  that  is  a  spherical 
portion  enclosed  by  the  arcs  of  great  circles. 

For  the  scenario  just  described  we  will  prove  a  lemma,  and  then  use  the  lemma 
to  produce  our  formula  for  u( 0, 0, 0,  T)  in  the  case  when  F  is  a  single  plane.  Let  S* 
be  the  angular  parameter  space  specifying  a  spherical  polygon  on  the  unit  sphere, 
and  let  B*(r,t)  C  B(r,t)  be  the  corresponding  spherical  polygon  within  the  sphere 
B{r,t).  The  boundary  dB*(r,t )  of  B*(r,t)  is  a  closed,  continuous,  and  piecewise 
smooth  curve  y(r,  t),  and  it  may  in  fact  be  a  single  great  circle.  In  any  case,  express¬ 
ing  the  boundary  y(r,  t)  as  a  union  Uj7j(r,  t)  of  smooth  curves,  we  use  dcr,;  =  rd(f>i 
to  represent  the  induced  Riemannian  measure  (differential  of  arc-length)  on  the 
component  7 i(r,t),  where  </>,;  is  an  angular  coordinate  along  the  component.  Fur¬ 
thermore,1  d/dxz  will  represent  the  Cartesian  direction  which  coincides  on  7 i(r,t) 
with  the  circle’s  outward  pointing  normal  as  a  component  of  dB*(r,t).  Along 
7 i(r,t)  the  vector  field  d/dx‘  points  perpendicularly  to  7 i(r,t)  and  also  tangent 
to  B(r,t).  Globally,  d/dxl  is  merely  the  normal  vector  field  for  some  foliation  of 
R3  into  R2  planes.  Let  Ad  represent  the  solid  past  null  cone  (or  conoid)  of  the 
spacetime  point  (0,  0,  0,T).  For  a  generic  time  t  <  T,  let  Mt  represent  the  closed 
portion  of  Ad  lying  to  the  future  of  time  level  t. 

Lemma:  Suppose  u  is  a  classical  solution  to  the  wave  equation  on  a  neighborhood 
of  A it..  Then,  with  $  the  solid  angle  subtended  by  B*(r,t),  we  have2 

$w(0,0,0,T)  [  ttt1 —  [  ^7 daldT 

i  Jt  T-t  Jli(T_TT)dxl 

(3.81)  [T  —  t)(ut)  +  9t{T  —  t)(u)  B,(T_t>ty 

The  result  can  be  shifted  to  a  generic  spatial  point  x  by  translation  invariance. 

To  prove  the  lemma,  we  first  note  that  the  equation  holds  in  the  t  — >  T~  limit. 
Indeed,  in  each  integral  over  "fi(T  —  r,  r)  the  apparently  singular  (T  —  r)_1  is 
canceled  by  a  (T  —  r)  in  the  dai  measure.  Therefore,  to  gather  the  result,  we  must 
simply  establish  that  the  right-hand  side  of  (3.81)  is  constant  in  t.  With  that  aim, 


1We  use  the  san  serif  x  to  indicate  that  the  d/dx‘  direction  need  not  be  one  of  the  fixed 
Cartesian  basis  directions:  d/dx,  d/dy,  d/dz. 

2For  w  =  w(x,  y,  z,t)  we  introduce  the  following  convention  for  (unnormalized)  angular  aver¬ 
ages: 

/,  w(r  sin  9  cos  (f),  r  sin  9  sin  </>,  r  cos  9 ,  t)dS, 

where  B*(r,t)  is  the  radius-r  spherical  portion  centered  at  the  origin  for  which  ( 9 ,  rf>)  E  <S*.  This 
average  does  not  use  the  proper  area  measure  r2dS  on  B*(r ,  t),  where  dS  =  sin  9d9dcj)  is  the  proper 
area  measure  on  the  unit— sphere.  By  choosing  not  to  incorporate  the  proper  area  measure  in  the 
definition  of  the  average,  we  ensure  that  r  — >  0"1”  limits  are  readily  taken. 
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we  consider  the  following  key  identity: 

(T-t)-1AS2u((T-t)v,t)  = 

(3-82)  ^  {(T  -  i)ut  ((T  -  t)u,  t)  +  A  [(T  -  t)u((T  -  t)u,  i)]  | , 

where  A52  is  the  S2  Laplacian  and  v  =  (sin^cos^,  sin0sin<(>,  cos0)  a  set  of  di¬ 
rection  cosines.  This  identity  is  nothing  more  than  the  wave  equation  itself,  here 
expressed  in  spherical  polar  coordinates.  Indeed,  note  that  the  left-hand  side  of 

(3.82)  is  symbolically 

(3.83)  ^  R(ut  +  uR  +  u/R)  =  R(utt  -  uRR  -  2u r/R). 

Since  the  angular  parameter  space  S*  does  not  depend  on  time,  we  may  integrate 
(3.82),  thereby  obtaining 

(T-t)-1  j  AS2U((T-t)v,t)dS  = 

(3.84)  IJ^T-  t)ut{(T  -  t)u,  t)  +  A  [(T  -  t)u{{T  -  t)u,  i)]  |  dS, 
or  in  our  more  compact  notation, 

(3.85)  {T-t)-\AS2u)B,(T_tt)  =  3t(T-f)<ut)B,(T_t^+at9T(T-f)<u)B,(T_M). 
By  Stokes’  Theorem,  the  term  on  left-hand  side  of  the  equation  integrates  to 

(3.86)  = 

that  is  precisely  minus  the  time  derivative  of  the  first  term  on  the  right-hand  side 
of  (3.81).  Whence  the  right-hand  side  of  (3.81)  is  indeed  constant  in  t,.  □ 

When  B*(T—t,t)  is  the  entire  sphere  B(T  —  t,t),  the  lemma  yields  the  standard 
spherical  means  formula 

u(0, 0,  0,  T )  = 

^JsM(T-t)v,t)dS  +  A  . 

We  see,  in  some  sense,  that  the  spherical  means  formula  holds  because  the  “bound¬ 
ary  of  a  boundary  is  zero”  [51].  See  [67]  for  another  derivation  of  (3.87).  We  remark 
that  the  lennna  above  can  also  be  established  via  a  generalization  of  Hadamard’s 
method  for  deriving  the  spherical  means  formula  [33,  18]. 

When  ry(T—t ,  t)  is  the  equatorial  great  circle  lying  in  the  plane  z  =  0,  the  lemma 
yields  a  hemispherical  means  formula , 

u(0, 0,  0,  T )  = 

1  fT  f27r 

—  —  /  /  uz((T  —  t)cos</>,  (T  —  r)sin<^,  Q,r)dcj)dT 

271"  Jt  Jo 

(3.88)  +P^j  Ut({T-t)v,t)dS+-^  j  u((T-t)u,t)dS  , 

with  S+  =  {(0,  <f>)  :  0  <  6  <  7t/2,  0  <  (f>  <  27r}  representing  the  angular  parameter 
space  specifying  the  northern  hemisphere.  If  we  take  t  =  0  as  the  initial  time  and 
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further  assume  that  the  initial  data  vanishes  for  z  >  0  (or  just  on  a  neighborhood 
of  B+(T,  0)),  then  the  last  formula  becomes 

1  fT  r2* 

(3.89)  u(0,  0, 0,  T)  =  —  —  /  /  uz((T  —  r)  cos  <f>,  (T  —  r)  sin  cf>,  0,  r)d(j>dT. 

2t  Jo  Jo 

A  similar  formula  will  always  follow  from  the  lemma,  provided  that  the  initial  data 
is  appropriately  chosen.  Since  the  2-derivative  uz  of  the  solution  also  obeys  the 
wave  equation,  the  last  equation  holds  with  u  replaced  by  uz  (again  subject  to  our 
assumption  about  initial  data).  Under  the  integral  sign,  one  can  then  exploit  the 
wave  equation  in  cylindrical  coordinates  to  derive  the  plane  boundary  condition 
described  in  [35]  as  expressed  in  terms  of  the  nonlocal  operator  in  (2.20),  although 
now  for  the  plane  2  =  0  rather  than  x  =  0. 

3.3.  Fast  evaluation  of  the  retarded  potential:  the  multilevel  plane  wave 
fast  time  domain  algorithm.  The  geometrical  flexibility  of  the  retarded  poten¬ 
tial  formulations  of  exact  boundary  conditions  is  obviously  attractive.  In  particular, 
unlike  the  formulations  we  have  presented  on  symmetric  boundaries,  they  allow  one 
to  use  a  computational  domain  of  minimal  size.  Let  us  then  consider  the  direct  dis¬ 
cretization  of  (3.78),  as  in  [25],  or  of  Teng’s  single-boundary  reformulation  (3.79) 
or  (3.80).  For  each  point  on  the  boundary  we  must  compute  an  integral  over  the 
boundary  of  data  extending  into  the  past.  By  our  scaling  assumptions  this  is  an 
0(1)  time  history  independent  of  T.  Thus  the  total  cost  of  a  direct  algorithm  is: 

(3.90)  Work  oc  A _5T,  Storage  oc  A~3. 

Thus,  although  the  storage  costs  are  comparable  to  those  required  by  the  volume 
solver,  the  work  requirements  are  excessive.  These  follow  from  the  dense  matrix 
multiplication  inherent  in  the  discretization  of  the  integrals.  The  analogous  problem 
arises  in  the  solution  of  frequency-domain  integral  equations  of  scattering  theory. 
For  the  frequency-domain  problem  there  has  been  an  intense  interest  in  the  inven¬ 
tion  of  fast  algorithms  to  compute  these  dense-matrix  multiplications.  Important 
examples  include  the  frequency-domain  fast  multipole  algorithm  (see  [17]  and  ref¬ 
erences  therein  for  the  mathematical  description  and  [15]  for  the  description  of  a 
large-scale  electromagnetic  scattering  code  which  uses  it)  and  the  equivalent  source 
method  (see  [12]  and  references  therein).  It  is  natural  to  attempt  to  apply  these 
methods  in  the  time  domain,  essentially  by  inverting  the  Fourier  transform.  The 
time-domain  version  of  the  fast-multipole  method  is  currently  the  best  developed 
algorithm  of  this  type  and  so  we  will  outline  it  below.  Equivalent  source  methods 
are  discussed  in  the  next  section. 

An  algorithm  for  evaluating  the  retarded  potential  based  on  fast-multipole  in¬ 
spired  ideas  is  the  multilevel  plane- wave  fast  time  domain  algorithm  (PWFTD)  of 
Michielssen  et  al  [23,  62].  The  details  of  its  implementation  are  somewhat  com¬ 
plex,  so  we  will  content  ourselves  with  an  overview,  referring  the  reader  to  the 
original  papers  for  more  details.  The  fundamental  ingredient  in  this  algorithm  is 
the  efficient  evaluation  of  space-time  localized  pieces  of  the  retarded  potential  in¬ 
tegral.  Consider  the  restriction  of  the  retarded  potential  integral  to  a  small  region 
of  space-time,  S  x  ( ts,tf ): 
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where  r  =  |x  —  x'|  is  the  distance  to  the  target  point,  x'.  Clearly,  this  contribution 
is  nonzero  only  in  the  time  interval  ( Ts,Tf )  with: 

(3.92)  Ts  =  ts  +  minr/c,  TV  =  1 1  +  maxr/c. 

xes  xes 

This  remarkable  property  of  solutions  of  the  wave  equation  and  Maxwell’s  equations 
in  three  space  dimensions  is  referred  to  as  the  strong  Huygens’  prinicple  or  the 
presence  of  lacunae. 

The  basis  of  the  PWFTD  is  the  representation  of  (3.91)  using  propagating  plane 
waves.3  The  mathematical  basis  for  the  time-domain  plane  wave  representations 
is  found  in  the  work  of  Heyman  [44],  He  points  out  the  fundamental  fact  that 
propagating  plane  wave  representations  are  not  causal;  as  soon  as  a  plane  wave  is 
“turned  on”  the  signal  is  felt  at  points  arbitrarily  remote  from  S.  Thus  the  true 
signal  requires  in  addition  evanescent  modes  which  precisely  cancel  these  noncausal 
signals.  However,  a  remarkable  result  of  [44]  is  that,  for  the  compactly  supported 
signals  considered  here,  regions  of  space-time  can  be  identified  where  only  the  prop¬ 
agating  waves  are  needed.  The  PWFTD  algorithm  employs  only  propagating  plane 
waves  to  evaluate  (3.91)  at  remote  locations  where  Heyman’s  analysis  shows  they 
are  sufficient.  The  outline  of  the  basic  two-level  algorithm  is  then  as  follows: 

i:  Expand  the  local  signal  into  a  discrete  set  of  plane  waves  with  directions 
appropriately  sampled  on  the  unit  sphere, 
ii:  Translate  the  planes  waves  to  remote  locations,  F. 
iii:  Evaluate  the  plane  waves  at  remote  nodes  within  F. 

This  basic  process  is  then  embedded  in  a  multilevel  framework.  The  final  result  is 
an  algorithm  requiring: 

(3.93)  WorkpwFTD  oc  A  3Tln  —  •  In  — , 

6  A 

which  is  formally  negligible  in  comparison  with  the  volume  solver.  However,  it  seems 
that  the  constants  are  larger  than  for  the  other  fast  methods  we  have  discussed. 

4.  Methods  Based  on  Equivalent  Sources 

Our  final  example  of  exact,  nonlocal  conditions  is  based  on  the  fact  that  solu¬ 
tions  to  Maxwell’s  equations  in  the  neighborhood  of  our  artificial  boundary  can  be 
represented  as  the  solution  of  the  forced  Maxwell  system  in  free  space,  with  sources 
distributed  in  the  region  between  any  scatterers  or  other  inhomogeneities  and  the 
artificial  boundary.  (For  example,  distributed  near  the  inner  surface  T  in  Figure 
1.)  An  algorithm  then  follows  from: 
i:  Finding  the  sources; 

ii:  Efficiently  evaluating  the  source  solution  at  the  boundary,  making  use  of 
the  strong  Huygen’s  principle  (the  existence  of  lacunae) . 

We  remark  that  the  retarded  potential  equation,  particularly  in  the  separated 
boundary  form  used  in  [25,  42],  can  be  viewed  as  a  special  case.  The  algorithms 
discussed  here  will  exploit  the  possibility  of  more  flexible  representations  to  derive 
efficient  algorithms. 


3A  curious  fact  about  the  algorithms  mentioned  here  is  that  they  employ  plane  wave  rather 
than  multipole  solution  representations  to  achieve  efficiency,  but  retain  the  word  multipole  in  their 
nomenclature  for  historical  reasons. 
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The  most  well-developed  equivalent  source  method  to  date  was  proposed  by 
Ryaben’kii  et  al  [58],  and  applied  to  Maxwell’s  equations  by  Tsynkov  in  [71].  We 
first  consider  their  construction  of  the  sources.  Define  a  cutoff  function, 
which  equals  one  near  and  beyond  the  artificial  boundary  and  zero  in  the  region 
containing  inhomogeneities;  /x  is  nonconstant  only  in  a  transition  region  which  ulti¬ 
mately  will  contain  the  sources.  In  applications  to  the  scalar  wave  equation  one  can 
simply  define  auxiliary  fields  by  multiplying  the  physical  fields  by  /x.  For  applica¬ 
tions  to  Maxwell’s  equations  this  direct  approach  is  problematic  as  the  subsequent 
time-decomposed  sources  fail  to  satisfy  the  continuity  relations.  Therefore,  a  some¬ 
what  more  involved  construction  is  advocated  in  [71]  which  we  will  only  outline 
here.  It  entails  the  construction  of  fields  W  and  V  satisfying: 

(4.94)  V  x  V  x  lb  =  B,  VxVxV  =  B, 
for  x  at  and  beyond  T'.  Then  set: 

(4.95)  B  =  VxVx  (fiW),  £  =  VxVx  (p,V). 


By  the  simple  application  of  the  product  rule  these  fields  satisfy  the  free  space 
Maxwell  system: 


(4.96) 

(4.97) 

as  well  as 

(4.98) 


dE 

— —  cVxB  =  -47 rj, 
at 

dB 

—  +  cV  x  E  =  -47rjm, 
at 


\7  ■  E  =  V  ■  B  =  0, 


(4.99)  V  •  j  =  V  •  jm  =  0. 

where  the  artificial  currents  are  defined  via  the  derivatives  of  /j  and  thus  are  nonzero 
only  in  the  transition  region.  The  purpose  of  the  indirect  construction  of  E  and  B 
is  to  guarantee  (4.98)-(4.99).  In  [71]  it  is  shown  that  W  and  V  may  be  determined 
only  from  the  knowledge  of  the  solution  on  T',  but  their  actual  construction  is  only 
described  in  a  special  case.  Thus  at  present  the  optimization  of  this  aspect  of  the 
algorithm  is  an  open  issue,  and  alternative  methods  are  being  studied  [72].  We 
emphasize  that  for  the  scalar  wave  equation  there  is  no  issue  here  as  the  auxiliary 
fields  can  be  defined  simply  through  multiplication  by  /x.  Of  course  practical  ap¬ 
plications  of  the  method  require  specific  choices  for  the  cutoff  function;  refer  to  the 
original  papers  [58,  71]  for  specific  examples. 

We  now  address  the  second  problem,  namely  the  efficient  evaluation  of  the  auxil¬ 
iary  fields  E  and  B  at  the  artificial  boundary  T'.  Recalling  that  these  fields  coincide 
with  E  and  B  on  T'  they  can  be  used  to  provide  exact  boundary  data  of  any  conve¬ 
nient  type.  Here  a  memory-efficient  algorithm  is  based  on  the  presence  of  lacunae. 
Consider  a  current  source  supported  in  the  interval  ( ts,tf ).  By  the  volume  Kirch- 
hoff  integral  we  have,  following  the  same  reasoning  as  in  the  preceding  section,  that 
the  solution  is  nonzero  at  T'  only  in  the  time  interval  U x'(ts,t/)  given  by  (3.92) 
where  x'  varies  over  I  '.  The  details  of  turning  this  fact  into  an  efficient  algorithm 
are  described  in  [59] .  A  sufficiently  smooth  partition  of  unity  in  time  is  introduced 
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to  express  the  sources  as  a  sum  of  terms  with  compact  time  support: 

(4.100)  j  =  'y  ] jk )  jm  —  ^  ]  jk,rn : 

k  k 

with  the  fcth  term  supported  in  (: )•  Define  Ek,  Bk  to  be  the  solution  of 
(4.96)-(4.97)  driven  by  the  fcth  source.  (Note  that  (4.99)  is  preserved  by  the  time 
partition.)  It  can  be  solved  on  the  finite  time  interval  (TSik,Tftk),  on  an  enlarged 
domain  with  simple  boundary  conditions.  In  [58]  periodic  boundary  conditions  are 
recommended.  Use  of  these  enable  the  use  of  Fourier  spectral  discretizations  which 
may  be  quite  efficient. 

Clearly,  the  computational  complexity  of  this  algorithm,  at  least  for  the  scalar 
wave  equation  where  construction  of  W  and  V  are  unnecessary,  will  be  of  the  same 
order  as  the  interior  solver,  (2.54).  As  such  it  may  not  be  as  efficient  as  some  of 
the  more  elaborate  constructions  discussed  earlier.  However,  again  at  least  when 
applied  to  the  scalar  wave  equation,  it  is  by  far  the  simplest  exact  method  to 
implement. 

More  recently,  Bruno  and  Hoch  [13]  have  developed  an  alternative  equivalent 
source  algorithm  for  the  scalar  wave  equation  which  has  the  potential  for  greater 
efficiency.  It  is  essentially  a  time-domain  version  of  Bruno  and  coworkers  fast 
algorithm  for  frequency  domain  scattering  (e.g.  [12]  and  references  therein).  Its 
essential  feature  in  comparison  with  the  method  described  above  is  the  use  of  simple 
sources  (monopoles  and  dipoles)  on  a  sparse,  regular  grid.  It  is  the  sparsity  of  the 
source  distribution  combined  with  the  use  of  FFTs  which  leads  to  the  potential 
savings.  The  computation  of  the  source  strengths  follows  from  the  least  squares 
approximation  to  the  solution  data  near  the  boundary. 

5.  Local  approximations 

Lastly  we  discuss  what  are  undobtedtly  the  most  often  used  techniques;  the 
perfectly  matched  layer  (PML)  and  local  radiation  boundary  conditions.  These 
methods  provide  geometric  flexibility  and  greater  potential  for  generalizations  to 
inhomogeneous  or  even  nonlinear  systems.  The  PML  in  particular  is  extremely 
simple  to  implement.  Although  these  methods  are  not  directly  based  on  exact 
formulations,  they  are  convergent,  albeit  often  nonuniformly  in  time.  As  such 
they  are  a  viable  alternative  to  the  methods  we  have  discussed,  particularly  if  the 
accuracy  requirements  are  not  too  stringent  and  the  solution  time,  T,  is  not  too 
large. 

We  begin  with  a  discussion  of  the  details  of  these  two  approaches,  as  well  as  an 
interesting  alternative  formulation  which  in  some  sense  unifies  them.  We  will  then 
discuss  their  convergence  properties  and  compare  their  complexity  to  the  algorithms 
implementing  nonlocal  formulations. 

5.1.  The  perfectly  matched  layer.  The  perfectly  matched  layer,  introduced  for 
Maxwell’s  equations  by  Berenger  [8],  is  an  absorbing  layer  with  a  reflectionless  inter¬ 
face  with  the  computational  domain.  Berenger’s  original  formulation  had  the  defect 
of  being  only  weakly  well-posed,  and  his  construction  was  somewhat  unintuitive. 
Subsequently,  a  clearer  understanding  of  PML  as  a  complex  coordinate  stretch¬ 
ing  emerged  [16].  Mathematically,  the  clearest  formulation  of  PMLs  for  Maxwell’s 
equations  has  been  given  by  Petropoulos  [54],  which  we  follow  here.  For  simplicity 
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we  display  his  model  for  a  Cartesian  layer  in  the  x  coordinate  direction.  In  [54] 
spherical  and  cylindrical  layers  are  also  developed. 

The  layer  equations  are  most  easily  discussed  in  the  frequency  domain,  where 
we  will  use  Laplace  rather  than  Fourier  transformations.  We  assume  the  layer  is 
located  in  x  £  (0,  L).  The  equations  then  follow  from  the  complex  coordinate 
stretching: 

(5101)  i  =  [  1  (’  +  7Tz)  iv’ 

where  s  is  the  dual  variable  to  time.  Maxwell’s  equations  then  become: 

(5.102)  SE-cVxB  =  0, 

(5.103)  SB  +  cVxE  =  0, 


where  V  x  is  obtained  by  replacing  the  x  derivatives  in  V  x  by: 


(5.104) 


1  s  +  a  d 
77  s  +  a  +  <j(x)  dx 


Time-domain  realizations  of  these  equations  are  obtained  by  viewing  the  trans¬ 
formed  system  as  an  anisotropic  dielectric  material.  The  layer  equations  then  are: 

dD 

(5.105)  — —  cV  x  H  =0, 

f)  Ft 

(5.106)  — +  cVx£=0, 


with  the  constitutive  relations: 

(5.107)  +  aEx  =  77  +  (a  +  cr)Dx^j  , 

(5.108)  ?7  +  («  +  cr)£’tan^  =  +  “A™ 

(5.109)  +  aHx  =  77  (^7^-  +  (a  +  , 

(5.110)  77  +  («  +  cr)17tan^  =  dg*an  +  ccBta„ 


Thus  to  implement  the  PML  one  only  needs  to  solve  an  additional  set  of  ordinary 
differential  equations.  Of  course  one  must  also  choose  the  parameters.  The  stretch¬ 
ing  parameter,  77  >  1,  is  often  omitted;  that  is  one  takes  77  =  1.  The  parameter 
a  >  0  is  called  the  complex  frequency  shift.  Often  it  is  set  to  zero,  but  choosing 
it  nonzero  yields  enhanced  long-time  stability  [7].  The  function  <r( x)  >  0  is  the 
absorption  parameter.  The  error  estimates  described  later  show  that  77  f0  a(p)dp 
controls  the  error.  Typically  a  =  <Joxq  is  used  with  q  chosen  so  that  the  fields  have 
sufficient  differentiability  properties  to  allow  differencing  across  the  layer  interface. 
However,  if  very  high  order  methods  are  used  this  may  be  a  restriction.  An  alterna¬ 
tive  is  to  develop  multidomain  formulations  with  characteristic  matching  across  the 
interface.  Then  one  can  choose  a  to  be  constant,  or  to  vanish  only  to  first  order. 
See,  for  example,  [22], 
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Geometric  flexibility  arises  from  the  implementation  of  the  PML  in  domains 
bounded  by  planar  faces.  Then,  in  addition  to  the  single-variable  layer  discussed 
above,  corner  layers  are  required.  These  are  derived  by  applying  the  coordinate 
stretching  in  all  three  variables.  See  the  appendix  of  [54]. 

Concerning  the  mathematical  properties  of  the  layers,  strong  well-posedness  and 
stability  can  be  established  [54,  4,  6].  We  note  that  in  [4]  more  general  PML 
formulations  are  derived  based  on  a  different  viewpoint.  These  generalizations  are 
necessary  for  the  treatment  of  anisotropic  materials. 


5.2.  Convergent  local  boundary  condition  sequences.  Finally  we  reconsider 
the  oldest  class  of  domain  truncation  methods,  local  radiation  boundary  condi¬ 
tion  sequences.  For  the  scalar  wave  equation  such  sequences  were  formulated  two 
decades  ago  by  Higdon  [45,  46].  They  have  been  revitalized  by  a  number  of  new 
developments  which  we  will  discuss  below.  These  are: 

i:  Development  of  new  auxiliary  variable  formulations  allowing  straightfor¬ 
ward  implementations  to  arbitrary  order; 
ii:  Construction  of  corner  compatibility  conditions  connecting  auxiliary  vari¬ 
ables  at  adjacent  faces  enabling  implementations  in  polygonal  domains; 
iii:  Adaptive  determination  of  boundary  condition  order; 
iv:  Proofs  of  spectral  convergence  with  increasing  order. 

Our  description  below  will  follow  [39].  We  note  that  parallel  developments  for 
the  scalar  wave  equation  are  reported  in  [40,  24,  37]. 

Consider  a  planar  artificial  boundary,  x  =  0.  Our  starting  point  is  a  representa¬ 
tion  of  the  solution  as  a  superposition  of  propagating  and  evanescent  plane  waves, 
derived  under  the  assumption  that  all  inhomogeneities  lie  in  the  left  half  plane 
x  <  — 6  for  some  6  >  0. 


i(x,y,z,t)  = 


4>  (ct  —  x  cos  9,y,z,  9)  d9 


(5.111) 


e  axT(y,z,t,a)da, 


where  u  is  any  Cartesian  held  component.  Following  the  treatment  of  analogous 
expressions  in  deriving  the  translation  formulas  for  the  fast  multipole  method  [26] , 
approximate  (5.111)  by  some  quadrature  rule: 

nP~  1  ne 

(5.112)  u(x,y,  z,t)  «  ^  Wj<&j(ct  —  x  cos  9j,y ,  z)  +  dje~<TjXT j(t,  y ,  z). 

j= o  j= i 

Local  boundary  conditions  with  np  +  ne  auxiliary  functions  can  now  be  constructed 
which  are  exact  on  this  approximate  representation  independent  of  the  unknown 
functions  and  Y,,-.  Recursively  define  for  j  =  0,  ...,np  —  1,  again  for  each 
Cartesian  component: 


(5.113) 


d 


d 


cos  9j  Q-t+CQ^)^=[  cos  -  CAT  )  Vh+i 


d 


d 


dt  dx 


with  tpo  =  u  and  for  j  =  1 ,ne: 

(5-114)  ^np+j-l  =  (°j  -  VV+J 


(5.115) 


^»P+ne  —  0. 
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Upon  proving  by  induction  that 
equations,  one  derives  evolution 
normal  characteristic  variables: 

(5.116)  Rj  =  E2lj 

•j  =  E3  ,j  ' 


'  B3 ,3 


>  v3 


We 


(5.117)  Sj  =  E3j  -  Wj  =  E3J  +  B2J. 

then  have  for  j  =  0, . . . ,  np  —  1: 

(1  +  cos6g^|±i  +  (1  -  cos  Oj)^- 

9Bij  dBiJ+1  dE\  j  dEij+i_ \ 

dz  dz  dy  dy  )' 


(5.118) 


dSj 


(5.119) 


(5.120) 


(5.121) 
and  for  j  - 

(5.122) 


(1  +  COS^O^i  +  (1  -  COS  Qj  )  ~7^r 

dBij  dBij+i  dEij  dEij+i\ 

dy  dy  dz  dz  J 

9vj+ 1 


dz 

dVj  „ 

I— +  (l-cos 


(1  +  cos^)-^ 

( dB\j  dBij+i  dE\j 
\  dz  dz 


■ft  \uvi+ 

'ej)~dT 

dEij+i  \ 


dy  dy 


n^W3  n  ^dw3  + 1 

COs9^~df  +  (1  “  cos^)  gt 
-  xv~  .  dE1J+1\ 


(1 

( ,  dBiJ+i  ,  ®Eij 
\  dy  dy  dz 


dz 


-1,..., ne: 


dRnp+j  dRHp+j- 1 

—QT  +  —di - 


Rrip+j) 


(dBu. 
k  dz 


dt, 

dBi,np+j  dEitnp+j-i  dE\ripjrj 

dz  dy  dy 


+  ■ 


dt  6 

,  x  f  dB in„+7  — 1  dB\  n 

(5.123)  — c  I - \  +  + - zrE±l 

V  oy  dv 


dBnp+j  9Snp+j-i 

~ *"  ^  C<JjWnp+j- 

dEi^np+j—i  dEiUp^.j 


Bnp+j) 


dz 


dz 


dy 

dVrip+j- 1  ,  dVrip+3  ,  /"!/■ 

- +  ^  +  co-j  (Kp+j-l 

^-£'l,np+j 

dy 


dt 

(5.124)  c  ("^-"p+J-1  +  ' 

\  dz  dz 


dt 

dEi.np+j_i 

dy 


Vnp+j) 


dwnp+i- 1  awno+7- 

J  +  ^-(Wn^-r  -  ^P+,) 


<9t 

Ir.  ioc'i  f®Bitnp+j-i  i  dBit.np+j  ^  dEitTlp+j_i  ^  dEi^Hp. 

(5-125)  ~CV  d~y  +  dy  +  ^ 
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In  addition  for  all  j : 
(5.126) 


dBld  (dE2, j  0E3j\ 
dt  °  V  dz  dy  )  ’ 


(5.127) 


dEhj  c(9B3J  dB2,J\ 
dt  \  dy  dz  )  ' 


We  see  that  the  structure  of  the  recursions  corresponds  to  the  directions  of  the 
characteristics.  For  the  outgoing  characteristics,  Rj  and  Sj ,  one  can  solve  for  the 
time  derivatives  with  increasing  j,  naturally  starting  with  the  time  derivatives  of 
Rq  and  Sq  which  can  be  computed  from  the  interior.  The  normal  fields  satisfy 
equations  which  are  uncoupled  in  j  and  which  thus  can  be  solved  individually.  For 
the  incoming  characteristics,  Vj  and  Wj,  on  the  other  hand,  one  can  solve  with 
decreasing  j.  We  then  determine  the  boundary  condition  by  setting: 


(5.128) 


Vrip+Ue  ~  Wrtp+ne  —  0. 


The  combination  (5.118)-(5.128)  thus  provides  a  recipe  for  computing  the  time 
derivatives  of  the  incoming  characteristic  variables  given  the  time  derivatives  of  the 
outgoing  variables.  Comparing  with  (2.12)-(2.13)  and  supposing  (as  we  always  have 
in  our  numerical  experiments)  that  9q  =  0  we  see  that  the  nonlocal  terms,  1Z Vo  and 
1Z Wq,  are  approximated  by: 


(5.129) 


RVo 


dEi:  i 
dy 


dBhl 

dz 


dE, 


dBi 


(5. no)  TZWo  gz  .  dy  . 

Despite  the  lengthy  description,  the  implementation  of  these  conditions  is  straight¬ 
forward.  In  fact  (5.118)-(5.127)  is  simply  a  hyperbolic  system  on  the  boundary 
which  can  be  discretized  using  whatever  scheme  is  used  in  the  interior.  In  [39] 
arbitrary-order  implementations  using  a  high-order  discontinuous  Galerkin  method 
are  demonstrated. 

To  use  these  conditions  on  polygonal  domains,  corner  and  edge  compatibility 
conditions  must  be  derived  to  provide  boundary  conditions  for  the  auxiliary  hyper¬ 
bolic  systems.  This  is  accomplished  in  [39]  in  two  space  dimensions  for  sequences 
with  ne  =  0.  The  construction  there  is  somewhat  ad  hoc;  it  depends  on  formally  in¬ 
troducing  doubly  indexed  auxiliary  variables  satsfying  the  recursions  on  both  faces, 
writing  down  the  large  system  of  equations  which  govern  them,  and  algebraically 
eliminating  all  space  derivatives.  Experiments  show  that  this  procedure  is  stable 
and  accurate  up  to  very  high  order  (over  100).  However,  it  is  not  yet  justified 
mathematically,  and  a  simpler  approach  would  be  desirable. 


5.3.  Implementations  via  optimal  grids:  a  link  between  PML  and  local 
boundary  condition  sequences.  Lastly  we  mention  an  interesting  connection 
between  PMLs  and  high-order  local  boundary  condition  sequences  for  the  scalar 
wave  equation  developed  by  Asvadurov  et  al.  [5]  and  used  later  in  [32].  The 
essential  idea  is  to  study  the  effective  discrete  Dirichlet-to-Neumann  maps  produced 
by  discretization  of  the  layer  equations.  One  then  realizes  that  for  a  fixed  finite 
difference  or  finite  element  discretization,  one  can  use  the  complete  freedom  of 
mesh  location,  including  the  possibility  of  choosing  a  complex  mesh,  to  control  the 
properties  of  this  map.  If  the  grid  is  chosen  to  agree  with  the  complex  grid  stretching 
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of  a  PML,  then,  of  course,  a  discretization  of  that  PML  is  produced.  However, 
other  choices  are  shown  to  correspond  to  approximations  akin  to  those  presented 
above.  In  [5]  a  particular  set  of  Oj’s  is  produced  corresponding  to  optimal  rational 
approximations  under  the  assumption  of  certain  distributions  of  propagating  waves. 
In  [32]  a  simpler,  nonstaggered,  grid  is  considered. 

We  will  show  how  the  ideas  presented  in  [32]  apply  to  the  recursive  formulation 
(5.113)-(5.115).  Applying  a  Laplace  transform  in  time,  the  basic  idea  is  the  recog¬ 
nition  that  if  we  treat  the  indices  in  these  equations  as  discrete  x  node  indices  then 
the  transformed  recursion  equations  can  be  rearranged  so  that  they  take  the  form: 


(5.131) 


(5.132) 


+ 

i  | 

(  di>j+i 

dipj 

f  2c  \ 

“  2  1 

i  dx 

+  ~d^ 

V  s  cos  0j  1 

i>j+ 1  -  i>j 

1  | 

\ 

(  9i>j+i 

di>j 

(  —  \ 

\ao  ) 

~  2  ' 

^  dx 

Thus  they  are  formally  equivalent  to  a  discretization  (via  the  box  scheme)  of  the 
identity  with  grid  spacings  and  This  establishes  a 

connection  with  PML  under  this  particular  discretization  as  one  would  simply  use 
different,  s-dependent  grid  spacings. 

Guddati  and  Lim  [32]  go  on  to  use  this  formal  relationship  to  very  simply  derive 
corner  compatibility  conditions;  one  simply  proceeds  as  with  PML  and  solves  using 
the  tensor-product  mapped  grid.  Clearly,  much  needs  to  be  done  to  establish  the 
mathematical  validity  of  this  rather  formalistic  construction.  It  does,  nonetheless, 
raise  interesting  issues  concerning  the  relationship  of  these  two  local  approaches 
after  discretization. 


5.4.  Accuracy  of  the  local  approximations.  A  direct  approach  to  assessing 
the  accuracy  of  the  local  approximations  is  to  return  to  the  Laplace  domain  and 
compute  the  reflection  coefficient.  We  assume,  as  in  the  derivation  of  (5.111), 
that  all  inhomogeneities  are  located  to  the  left  of  x  —  —5  and  that  the  artificial 
boundary  is  the  plane  x  =  0.  Estimating  the  error  in  terms  of  the  reflection 
coefficient  is  a  straightforward  application  of  Parseval’s  relation;  see  [35,  36,  39]  for 
details.  Precisely,  for  fixed  tangential  wave  numbers  k  they  are  given  in  terms  of: 


(5.133) 


max  \R(s,k)\. 
St  «=T-! 


For  the  PML  reflection  coefficients  are  computed  in  many  places,  though  often 
only  displayed  in  the  propagating  mode  regime.  Here  we  take  r)  =  1,  suppose  a  layer 
of  width  L  terminated  by  normal  characteristic  boundary  conditions,  and  define: 


(5.134) 

After  some  straightforward  algebra  we  find: 


a  =  L  1  <j{t)(It. 

Jo 


( (s2  +  \k\2)1/2  -  s\ 

V(s2  +  |fc|2)1/2  +  sy 


(5.135) 


-ffpML 


e 


2L(S2  +  |fc|2)1/2(l+lff^) 
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Similarly,  one  can  compute  the  reflection  coefficient  for  the  local  boundary  con¬ 
dition  sequences  [39] : 


(5.136) 

where 


f?BC  —  Np  ■  Kjj 


/  (s2  +  |fc| 2)1/2  -  s\ 

Us2  +  |fc|  2)1/2+iJ’ 


(5.137) 


=  fnfj1  cos  6 jS  (s2  +  |fc|2)1/2\ 

l  cos  +  (^2  +  Ifcp)1/2  y 


(5.138) 


he 


f^^-(s2  +  \k\2)1/2\ 
\j=\  ^+(§2  +  m^j  ■ 


We  first  note  that  the  restriction  to  finite  times  is  necessary;  if  we  take  s  imag¬ 
inary  the  maximum  reflection  coefficients  are  one.  This  difficulty  exists  even  for 
the  nonlocal  conditions  on  planar  and  cylindrical  boundaries  and  is  responsible  for 
the  InT  terms  in  the  complexity  estimates.  Then  for  fixed  T  and  assuming  a  ban- 
dlimited  signal,  fc  <  A-1,  it  is  clear  that  both  methods  converge.  In  particular 
-ffpML  — »  0  as  Ld  — >  oo  since: 


(5.139) 


min  3? 

8ts=T-! 


(s2  +  Ifcl2)1/2^ 
s  +  a  1 


>  0. 


Similarly,  Rbc  — >  0  as  np  — >  oo  as  each  term  in  the  definition  of  up  is  strictly 
smaller  than  one. 

Rather  than  attempting  to  build  analytic  error  estimates  out  of  these  expressions 
we  will  simply  compute  parameter  values  needed  to  meet  various  tolerances.  We 
consider  three  cases:  T  =  10,  A  =  10-1;  T  =  20,  A  =  2  x  1CP2;  and  T  =  100, 
A  =  10-2.  In  the  case  of  the  PML  we  set  a  =  0.1,  a  =  1,  and  only  vary  L.  In  the 
plots  we  assume  that  the  number  of  points  in  the  layer  is  proportional  to  A_1L.  In 
practice  one  can  decrease  the  resolution  within  the  layer  so  that  we  are  somewhat 
overpredicting  the  number  of  points  required,  but  we  won’t  attempt  to  quantify 
this  effect. 

For  the  boundary  conditions  we  consider  two  choices. 

Pade  Parameters: 


(5.140)  cos  9 j  —  1 ,  j  =  0, . . . ,  rip  —  13  ne  =  0, 
Gauss-Radau-Rokhlin-Yarvin  (GRRY)  Parameters: 

(5.141)  dj  =  7r(Cj4+1),  j  =  0, . . .  ,np  —  1, 

where  Cj  are  the  left  endpoint  Gauss- Radau  nodes  on  [— 1, 1]  and 

(5.142)  aj=(3\dj,  j  =  Q,...,ne, 

where  dj  are  the  Yarvin-Rokhlin  nodes  [75].  (The  Yarvin-Rokhlin  nodes  are  tabu¬ 
lated  in  the  f77  subroutine  wts500.f  available  at  www.netlib.org/pdes/multipole/.) 
For  the  experiments  we  chose  np  =  ne  and  /3  =  5. 

In  Figure  2  we  plot  the  maximum  of  the  reflection  coefficient  as  a  function  of  the 
degrees  of  freedom  in  the  boundary  treatment  (terms  in  the  boundary  condition  or 
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points  in  the  layer).  These  can  be  approximately  fit  with  the  following  exponential 
convergence  models: 


(5.143) 


errpML 

errpade 

eri'GRRY 


oc 


oc 


oc 


Comparitive  numerical  experiments  are  presented  in  [39] .  There  we  see  that  the 
actual  errors  are  typically  an  order  of  magnitude  or  two  smaller  than  predicted  by 
the  maximum  reflection  coefficient.  We  also  show  that  the  PML  results  can  be 
improved  on  by  coarsening  the  the  layer  resolution.  Nonetheless,  the  results  do 
follow  the  trends  predicted  by  this  analysis;  the  long-time  error  is  worse  for  the 
traditional  local  Pade  boundary  conditions,  using  the  PML  results  in  significant 
improvements,  and  use  of  the  GRRY  parameters  leads  to  significant  improvements 
still. 

Accepting  the  estimates  in  (5.143)  we  can  finally  estimate  the  complexity  of  the 
local  boundary  treatments: 


_7  2  1 
WorkpML  oc  A  2T2ln-, 
e 

(5.144)  Workpade  oc  A”4T2  In  -, 

e 

WorkcRRY  oc  A_3Tln  -  •  In  (A-1T), 
e 


StoragePML  oc  \~iT^  In  A 
StoragePade  a  A”3Tln  - 
StorageGRRY  oc  A-2  In  A  •  In  (A_1T). 


We  see  that  the  PML  approach  is  acceptable  except  for  T  very  large.  The  use  of 
traditional  local  boundary  conditions  is  only  acceptable  for  T  =  0(1).  The  new 
local  boundary  conditions,  on  the  other  hand,  yield  complexity  estimates  compa¬ 
rable  to  the  nonlocal  boundary  conditions.  Thus  if  it  would  be  possible  to  extend 
the  construction  of  corner  compatibility  conditions  to  this  case  it  seems  they  would 
provide  a  fairly  complete  solution.  However,  at  present  the  corner  compatibility 
conditions  have  only  been  constructed  for  the  case  ne  =  0. 

Lastly  we  note  that  exact  reflection  formulas  have  recently  been  derived  by  Diaz 
and  Joly  [20,  21]  and  de  Hoop  et  al  [19].  They  study  the  scalar  wave  equation  and 
use  the  Cagniard-de  Hoop  method.  This  somewhat  restricts  the  parametrizations 
they  can  study.  In  particular  only  the  cases  ne  =  0  for  boundary  condition  se¬ 
quences  and  a  =  0  for  PML  are  treated.  In  these  cases  the  results  are  in  agreement 
with  those  stated  above. 

Our  numerical  experiments  for  time-dependent  problems  in  waveguides  confirm 
the  T-dependence  in  the  complexity  estimates  above.  However,  for  spherical  bound¬ 
aries  they  may  be  overly  pessimistic.  A  detailed  analysis  for  spherical  PMLs  in  the 
frequency  domain  has  been  given  by  Bao  and  Wu  [10]  However,  the  frequency- 
dependence  of  their  estimates  precludes  their  direct  application  to  time-domain 
problems.  Recently,  Chen  has  completed  a  time-domain  analysis  of  spherical  PMLs 
for  the  scalar  wave  equation  [14],  but  it  is  unclear  at  present  if  the  T-dependence 
in  his  estimates  is  optimal. 
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Reflection  coefficient  for  PML 


Reflection  coefficient  for  Pade  parameters 


Figure  2.  Maximum  reflection  coefficients  for  various  local  do¬ 
main  truncation  techniques. 
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6.  Epilogue:  towards  the  ultimate  solution 

We  have  seen  that  there  are  now  a  number  of  techniques  that  provide  completely 
satisfactory  solutions  to  the  time-domain  radiation  boundary  condition  problem  in 
electromagnetics  for  a  wide  range  of  situations.  At  present,  none  can  be  called 
optimal  for  all  cases  considered.  Although  the  nonlocal  conditions  are  extremely 
efficient,  they  lack  geometric  flexibility.  In  addition,  they  have  not  been  generalized 
to  treat  the  important  case  of  mutiple  media  extending  to  infinity,  though  it  is  likely 
that  some  extensions  in  this  direction  are  possible.  Use  of  the  retarded  potential  in 
conjunction  with  the  PWFTD  does  have  favorable  complexity  estimates,  but  prac¬ 
tical  experience  with  the  algrithm  shows  that  there  is  substantial  computational 
overhead.  The  perfectly  matched  layer  does  possess  geometric  flexibility  and  ex¬ 
tensibility  to  more  complex  models.  However,  its  accuracy  can  suffer  in  long  time 
simulations  and  the  issues  associated  with  optimal  numerical  implementations  are 
not  easy.  Local  conditions  based  on  the  GRRY  nodes  show  the  most  promise,  but 
their  implementation  in  domains  with  corners,  which  is  necessary  if  they  are  to  be 
made  geometrically  flexible,  has  yet  to  be  demonstrated. 

We  finish  by  pointing  out  an  interesting  mathematical  result  due  to  Warchall 
[73].  We  state  it  for  the  scalar  wave  equation  with  some  simplifying  assumptions; 
a  more  general  version  is  proven  in  [73]. 

Theorem  1  (Warchall).  Let  Ll  C  R”  be  an  open  convex  set.  Let  f(x,t),Uo(x),Vo(x) 
be  sufficiently  smooth  and  compactly  supported  in  SV  C  fl.  Finally,  let  u  satisfy 
□u  =  /,  u(x,  0)  =  uq{x),  Qfjr(x,0)  =  Vo(x).  Suppose  x  €  dfl  and  At  is  such  that 
cAt  <  dist(x,  $T).  Then  if  u(x,t)  and  Qjf(x,t)  vanish  for  all  x  €  fl  satisfying 
\x  —  x\  <  cAt,  we  may  conclude  that  u(x,  t  +  At)  =  0. 

The  direct  interpretation  of  Warchall’s  Theorem  is  as  follows.  If  we  choose  a 
convex  artificial  boundary  such  that  the  support  of  the  data  are  located  some  small 
distance  away  and  a  time  step  so  that  the  a  priori  domain  of  dependence  of  boundary 
points  over  a  time  step  doesn’t  intersect  with  the  support  of  the  data,  then  any  two 
solutions  produced  by  (possibly  different)  data  with  the  same  support  which  agree 
on  the  restriction  of  the  domain  of  dependence  to  the  computational  domain  will 
have  the  same  updated  values  on  the  boundary. 

We  would  like  to  translate  this  local  uniqueness  result  to  a  formula  showing  how 
to  use  the  data  inside  U  to  update  the  solution.  If  a  stable,  local,  exact  update 
formula  could  be  found,  it  would  clearly  represent  the  ultimate  solution  to  the 
radiation  boundary  condition  problem. 
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